Awesome test case for demonstrating the bug! The implementation doesn't make sense to me though.

w[end] += σ - sum(w) could end up being negative depending on the current value of w[end], which breaks the definition of simplex. Furthermore, it doesn't seem very rigorous to subtract from an arbitrary element in the vector.

It seems like this is a problem of insufficient precision, so I propose that we instead convert to BigFloat and back in the method.

Proposed implementation:

function  simplex(x::AbstractVector{T},  ::Real) where {T<:Real}
    _x = convert(Vector{BigFloat}, x)
    _sigma  = convert(BigFloat,  )
    n = length(_x)
    mu = sort(_x; rev=true)
    rho = maximum((1:n)[ -(cumsum(mu).-_sigma)./(1:n).>zero(BigFloat)])
    theta = (sum(mu[1:rho]) - _sigma) / rho
    _w = max.(_x .- theta, zero(BigFloat))
    w = convert(Vector{T}, _w)
    return w
end

Going to close this PR without merging as the proposed implementation passes @hopeyen's test and stays true to the definition of a simplex projection.